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1 Introduction 

In these lectures I cover a number of topics in cosmological data analysis. I concentrate on 
general techniques which are common in cosmology, or techniques which have been developed 
in a cosmological context. In fact they have very general applicability, for problems in which 
the data are interpreted in the context of a theoretical model, and thus lend themselves to 
a Bayesian treatment. 

We consider the general problem of estimating parameters from data, and consider how 
one can use Fisher matrices to analyse survey designs before any data are taken, to see 
whether the survey will actually do what is required. We outline numerical methods for esti- 
mating parameters from data, including Monte Carlo Markov Chains and the Hamiltonian 
Monte Carlo method. We also look at Model Selection, which covers various scenarios such 
as whether an extra parameter is preferred by the data, or answering wider questions such 
as which theoretical framework is favoured, using General Relativity and braneworld gravity 
as an example. These notes are not a literature review, so there are relatively few references. 
Some of the derivations follow the excellent notes of Licia Verde [23] and Andrew Hamilton 

After this introduction, the sections are: 

• Parameter Estimation 

• Fisher Matrix analysis 

• Numerical methods for Parameter Estimation 

• Model Selection 

1.1 Notations: 

• Data will be called x, or x, or Xi, and is written as a vector, even if it is a 2D image. 

• Model parameters will be called 0, or 9 or 9 a 

1.2 Inverse problems 

Most data analysis problems are inverse problems. You have a set of data x, and you wish 
to interpret the data in some way. Typical classifications are: 

• Hypothesis testing 

• Parameter estimation 



• Model selection 
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Cosmological examples of the first type include 

• Are CMB data consistent with the hypothesis that the initial fluctuations were gaussian, 
as predicted (more-or-less) by the simplest inflation theories? 

• Are large-scale structure observations consistent with the hypothesis that the Universe 
is spatially fiat? 

Cosmological examples of the second type include 

• In the Big Bang model, what is the value of the matter density parameter? 

• What is the value of the Hubble constant? 

Model selection can include slightly different types (but are mostly concerned with larger, 
often more qualitative questions): 

• Do cosmological data favour the Big Bang theory or the Steady State theory? 

• Is the gravity law General Relativity or higher-dimensional? 

• Is there evidence for a non-flat Universe? 

Note that the notion of a model can be a completely different paradigm (the first two 
examples), or basically the same model, but with a different parameter set. In the third 
example, we are comparing a flat Big Bang model with a non-flat one. The latter has 
an additional parameter, and is considered to be a different model. Similar problems exist 
elsewhere in astrophysics, such as how many absorption-line systems are required to account 
for an observed feature? 

These lectures will principally be concerned with questions of the last two types, param- 
eter estimation and model selection, but will also touch on experimental design and error 
forecasting. Hypothesis testing can be treated in a similar manner. 



2 Parameter estimation 

We collect some data, and wish to interpret them in terms of a model. A model is a theoretical 
framework which we assume is true. It will typically have some parameters in it, which 
you want to determine. The goal of parameter estimation is to provide estimates of the 
parameters, and their errors, or ideally the whole probability distribution of 9, given the 
data x. This is called the posterior probability distribution i.e. it is the probability that the 
parameter takes certain values, after doing the experiment (as well as assuming a lot of other 
things): 

p(0|x). (I) 

This is an example of RULE Start by thinking about what it is you want to know, and 
write it down mathematically. 

From p{6\x) one can calculate the expectation values of the parameters, and their errors. 
Note that we are immediately taking a Bayesian view of probability, as a 'em degree of 
belief, rather than a frequency of occurrence in a set of trials. 

2.1 Forward modelling 

Often, what may be easily calculable is not this, rather the opposite, p(x|0)[^] The opposite is 
sometimes referred to as forward modelling - i.e. if we know what the parameters are, we can 

1 There is no rule n: n > 1. 

2 If you are confused about p(A\B) and ^(51^4) consider if A=pregnant and B=female. p(A\B) is 
a few percent, p(B\A) is unity 
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compute the expected distribution of the data. Examples of forward modelling distributions 
include the common ones - binomial, Poisson, gaussian etc. or may be more complex, such 
as the predictions for the CMB power spectrum as a function of cosmological parameters. 
As a concrete example, consider a model which is a gaussian with mean /x and variance a 2 . 
The model has two parameters, — (/i, a), and the probability of a single variable x given 
the parameters is 

(*-m) 21 



P(x\0) = _ exp 

V27TCT 



2a 2 



(2) 



but this is not what we actually want. However, we can relate this to p(0\x) using Bayes' 
Theorem, here written for a more general data vector x: 

pi e\ x ) = p W= p WM°\ (3) 

p(x) p(x) 

• p(0\x) is the posterior probability for the parameters. 

• p(x\0) is called the Likelihood and given its own symbol £(x; 0). 

• p(6) is called the prior, and expresses what we know about the parameters prior to 
the experiment being done. This may be the result of previous experiments, or theory 
(e.g. some parameters, such as the age of the Universe, may have to be positive). In the 
absence of any previous information, the prior is often assumed to be a constant (a 'flat 
prior'). 

• p(x) is the evidence. 



For parameter estimation, the evidence simply acts to normalise the probabilities, 

p(x)= J d8p(x\8)p(0) (4) 

and the relative probabilities of the parameters do not depend on it, so it is often ignored 
and not even calculated. 

However, the evidence does play an important role in model selection, when more than 
one theoretical model is being considered, and one wants to choose which model is most 
likely, whatever the parameters are. We turn to this later. 

Actually all the probabilities above should be conditional probabilities, given any prior 
information I which we may have. For clarity, I have omitted these for now. / may be the 
result of previous experiments, or may be a theoretical prior, in the absence of any data. In 
such cases, it is common to adopt the principle of indifference and assume that all values 
of the parameter(s) is (are) equally likely, and take p(0)=constant (perhaps within some 
finite bounds, or if infinite bounds, set it to some arbitrary constant and work with an 
unnormalised prior). This is referred to as a flat prior. Other choices can be justified. 

Thus for flat priors, we have simply 

p(0|x) oc £(x;0). (5) 

Although we may have the full probability distribution for the parameters, often one simply 
uses the peak of the distribution as the estimate of the parameters. This is then a Maximum 
Likelihood estimate. Note that if the priors are not flat, the peak in the posterior p(0|x) is 
not necessarily the maximum likelihood estimate. 

A 'rule of thumb' is that if the priors are assigned theoretically, and they influence the 
result significantly, the data are usually not good enough. (If the priors come from previous 
experiment, the situation is different - we can be more certain that we really have some prior 
knowledge in this case). 

Finally, note that this method does not generally give a goodness-of-fit, only relative 
probabilities. It is still common to compute % 2 at this point to check the fit is sensible. 
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2.2 Updating the probability distribution for a parameter 

One will often see in the literature forecasts for a new survey, where it is assumed that we 
will know quite a lot about cosmological parameters from another experiment. Typically 
these days it is Planck, which is predicted to constrain many cosmological parameters very 
accurately. Often people 'include a Planck prior'. What does this mean, and is it justified? 
Essentially, what is assumed is that by the time of the survey, Planck will have happened, 
and we can combine results. We can do this in two ways: regard Planck+survey as new data, 
or regard the survey as the new data, but our prior information has been set by what we 
know from Planck. If Bayesian statistics makes sense, it should not matter which we choose. 
We show this now. 

If we obtain some more information, from a new experiment, then we can use Bayes' 
theorem to update our estimate of the probabilities associated with each parameter. The 
problem reduces to that of showing that adding the results of a new experiment to the 
probability of the parameters is the same as doing the two experiments first, and then 
seeing how they both affect the probability of the parameters. In other words it should not 
matter how we gain our information, the effect on the probability of the parameters should 
be the same. 

We start with Bayes' expression for the posterior probability of a parameter (or more 
generally of some hypothesis), where we put explicitly that all probabilities are conditional 
on some prior information /. 

Let say we do a new experiment with new data, x'. We have two ways to analyse the new 
data: 

• Interpretation 1: we regard x' as the dataset, and x/ (means x and I) as the new prior 
information. 

• Interpretation 2: we put all the data together, and call it x'x, and interpret it with the 
old prior information /. 

If Bayesian inference is to be consistent, it should not matter which we do. 

Let us start with interpretation 1. We rewrite Bayes' theorem, equation ^ by changing 
datasets x — > x', and letting the old data become part of the prior information / — > I' = x/. 
Bayes' theorem is now 

p(x'|xj) 

We now notice that the new prior in this expression is just the old posteriori probability 
from equation ([6]), and that the new likelihood is just 

™ = f^ 

Substituting this expression for the new likelihood: 

p(0|x/Mx'x|0I) 
P[ 1 ' - p(x'|x/)p(x|0/) ' W 

Using Bayes' theorem again on the first term on the top and the second on the bottom, we 
find 

p(x'|xi>(x|I) 

and simplifying the bottom gives finally 



(11) 
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which is Bayes' theorem in Interpretation 2. i.e. it has the same form as equation ([6]), the 
outcome from the initial experiment, but now with the data x replaced by x'x. In other 
words, we have shown that x — > x' and / — > xl is equivalent tox-> x'x. This shows us that 
it doesn't matter how we add in new information. Bayes' theorem gives us a natural way of 
improving our statistical inferences as our state of knowledge increases. 



2.3 Errors 



Let us assume we have a posterior probability distribution, which is single-peaked. Two 
common estimators (indicated by a hat: 0) of the parameters are the peak (most probable) 
values, or the mean, 



6 = J d66p(6\x). (12) 

An estimator is unbiased if its expectation value is the true value 9q: 

(6) = O . (13) 

Let us assume for now that the prior is flat, so the posterior is proportional to the likelihood. 
This can be relaxed. Close to the peak, a Taylor expansion of the log likelihood implies that 
locally it is a mutivariate gaussian in parameter space: 

In L(x; 0) = In L(x; 9 ) + \ (0 a - 8 0a ) (ty _ 0Qfj ) + . . . (14) 



or 

L(x; 9) — £(x; O ) ex P ' 



(15) 



o2 I r 

The Hessian matrix Y\ a p = — q 9 gg controls whether the estimates of 9 a and Op are corre- 
lated or not. If it is diagonal, the estimates are uncorrelated. Note that this is a statement 
about estimates of the quantities, not the quantities themselves, which may be entirely in- 
dependent, but if they have a similar effect on the data, their estimates may be correlated. 
Note that in cases of practical interest, the likelihood may not be well described by a mul- 
tivariate gaussian at levels which set the interesting credibility levels (e.g. 68%). We turn 
later to how to proceed in such cases. 



2.4 Conditional and marginal errors 



If we fix all the parameters except one, then the error is given by the curvature along a line 
through the likelihood (posterior, if prior is not flat): 

^"conditional, a y u (^^) 

V Hon 

This is called the conditional error, and is the minimum error bar attainable on 9 a if all the 
other parameters are known. It is rarely relevant and should almost never be quoted. 



2.5 Marginalising over a gaussian likelihood 



The marginal distribution of B\ is obtained by integrating over the other parameters: 

p(0i) = J d0 2 ...dO N p{0) (17) 

a process which is called marginalisation. Often one sees marginal distributions of all pa- 
rameters in pairs, as a way to present some complex results. In this case two variables are 
left out of the integration. 
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If you plot such error ellipses, you must say what contours you plot. If you say you 
plot lcr and 2cr contours, I don't know whether this is for the joint distribution (i.e. 68% 
of the probability lies within the inner contour), or whether 68% of the probability of a 
single parameter lies within the bounds projected onto a parameter axis. The latter is a lcr, 
single-parameter error contour (and corresponds to Ax 2 — 1), whereas the former is a lcr 
contour for the joint distribution, and corresponds to A\ 2 = 2.3. 

Note that Ax 2 = X 2 ~~ x 2 (minimum) , where 

X 2 = £ (18) 

i 0i 

for data Xi with /Ltj = (x^ and variance of. If the data are correlated, this generalises to 

\" »< l^ ( './ i ->-J I<j> (19) 

ij 

where = ((xj - m)(xj - fij)). 

For other dimensions, see Table 1, or read... 

The Numerical Recipes bible, chapter 15.6 [19j 

Read it. Then, when you need to plot some error contours, read it again. 

Table 1. A\ 2 for joint parameter estimation for 1, 2 and 3 parameters. 



a 


P 


M=l 


M=2 


M=3 


lcr 


68.3% 


1.00 


2.30 


3.53 


2cr 


95.4% 


4.00 


6.17 


8.02 


3cr 


99.73% 


9.00 


11.8 


14.2 



Note that some of the results I give assume the likelihood (or posterior) is well- 
approximated by a multivariate gaussian. This may not be so. If your posterior is a single 
peak, but is not well-approximated by a multivariate gaussian, label your contours with the 
enclosed probability. If the likelihood is complicated (e.g. multimodal), then you may have 
to plot it and leave it at that - reducing it to a maximum likelihood point and error matrix 
is not very helpful. Not that in this case, the mean of the posterior may be unhelpful - it 
may lie in a region of parameter space with a very small posterior. 

A multivariate gaussian likelihood is a common assumption, so it is useful to compute 
marginal errors for this rather general situation. The simple result is that the marginal error 
on parameter 9 a is 

o a = ^(H- 1 )^. (20) 

Note that we invert the Hessian matrix, and then take the square root of the diagonal 
components. Let us prove this important result. In practice it is often used to estimate 
errors for a future experiment, where we deal with the expectation value of the Hessian, 
called the Fisher Matrix. 

<21) 

We will have much more to say about Fisher matrices later. The expected error on 9 a is 
thus 

<r a = y/ (F _1 ) aQ ;. (22) 

It is always at least as large as the expected conditional error. Note: this result applies for 
gaussian-shaped likelihoods, and is useful for experimental design. For real data, you would 
do the marginalisation a different way - see later. 

To prove this, we will use characteristic functions. 



Characteristic functions 
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In probability theory the Fourier Transform of a probability distribution function is known as 
the characteristic function. For a multivariate distribution with N parameters, it is denned 
by 

ik0 



with reciprocal relation 



tf>(k) = [d N 6p(9)e- i] 

r r1 N k 



k-0 



(23) 



(24) 



(note the choice of where to put the factors of 2tt is not universal). Hence the characteristic 
function is also the expectation value of e : 



m = < e - ik -0>. 



(25) 



Part of the power of characteristic functions is the ease with which one can generate all 
of the moments of the distribution by differentiation: 



_d{-i\L a ) n ° ...d(-ik^f> 
This can be seen if one expands (f>(k) in a power series, using 



k=0 



exp(a) =22^1 ' 



(26) 



(27) 



giving 



0(k) = i - ik • {6) - k «M^ 



dp 



Hence for example we can compute the mean 

00(k) 



(0a) = 



d(-ik a ) 



k=0 



and the covariances, from 



a 9p) — 



d(-ik a )d(-ik[3) 



k=0 



(28) 



(29) 



(30) 



(Putting a = p yields the variance of a after subtracting the square of the mean) 



2.6 The expected marginal error on 6 a is y/ (F — 1 ) aa 

The likelihood is here assumed to be a multivariate gaussian, with expected hessian given 
by the Fisher matrix. Thus (suppressing ensemble averages) 



L(0) 



(27r) M /2VdcTF 



exp 



1 



e 1 fo 



(31) 



where T indicates transpose, and for simplicity I have assumed the parameters have zero 
mean (if not, just redefine as the difference between 6 and the mean). We proceed by 
diagonalising the quadratic, then computing the characteristic function, and compute the 



covariances using equation ( 30 1 . This is achieved in the standard way by rotating the pa- 
rameter axes: 

* = R0 (32) 

for a matrix R. Since F is real and symmetric, R is orthogonal, R" 1 = R T . Diagonalising 
gives 
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9 T F9 = * T RFR T *, 
and the diagonal matrix composed of the eigenvalues of F 

A = RFR T , 

Note that the eigenvalues of F are positive, as F must be positive-definite. 

The characteristic function is 
1 



fk) = 



(2 7 r) M /2 v / detF 



J d M $>exp (-^ T AF*j cxp(-ik T R T *) 



(33) 
(34) 

(35) 



where we exploit the fact that the rotation has unit Jacobian to change d M 9 to d M ^. If we 
define K = Rk, 



(36) 



and since A is diagonal, the first exponential is a sum of squares, which we can integrate 
separately, using 



dipexp(-Ai; 2 /2)exp(-iKi;) = y / 2ir/Acxp[-K 2 /{2A)}. 



(37) 



All multiplicative factors cancel (since the rotation preserves the eigenvalues, so det(F) 
n^la), and we obtain 



m = exp ^-Y^K 2 /{2A l )^j = exp (-^K T A^Kj = cxp 



-k T F- x k 



(38) 



where the last result follows from A K = k r (R y yl _i R)k = k^ F^k. 



Having obtained the characteristic function, the result (22 1 follows immediately from 



equation (30). 



2.7 Marginalising over 'amplitude' variables 



It is not uncommon to want to marginalise over a nuisance parameter which is a simple 
scaling variable. Examples include a calibration uncertainty, or perhaps an unknown gain in 
an electronic detector. This can be done analytically for a gaussian data space: 



L(0;x) 



1 



(2 7 r) 7V / 2 VdetC 



exp 



i^(x i -xf)Cr 1 (, 



(39) 



where we have N data points with covariance matrix Qj = ((x, — x*' l )(x :) — x'- )), and 
the model data values x* h and C depend on the parameters. Let us now assume that the 
covariance matrix does not depend on the parameters, so the parameter dependence is 
only through x th (6), and further we assume that one of the parameters simply scales the 
theoretical signal, i.e. x (9) = Ax th {0\A = 1). If we want the likelihood of all the other 
parameters (call this L'(9';x), with one fewer parameter), marginalised over the unknown 
A, then we can integrate: 



L(0';x) = J dAL(9-x)p(A) 



(2 7 r) Ar / 2 VdcTC 



dA exp 



-^(x i -Axf)Cr. 1 (x J --Axf) 



p(A) (40) 



where p(A) is the prior for A. If we take a uniform (unnormalised!) prior on A between 
limits ±oo, and the theoretical x are now at A = 1, then we can integrate the quadratic. 
You can do this. 



Exercise: do this. 
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3 Fisher Matrix Analysis 

This has been adapted from [H] (hereafter TTH). 

How accurately can we estimate model parameters from a given data set? This question 
was basically answered 60 years ago |5], and we will now summarize the results, which are 
both simple and useful. 

Suppose for definiteness that our data set consists of N real numbers x%, X2,-..,xn, 
which we arrange in an iV-dimensional vector x. These numbers could for instance denote 
the measured temperatures in the N pixels of a CMB sky map, the counts-in-cells of a galaxy 
redshift survey, TV coefficients of a Fourier expansion of an observed galaxy density field, 
or the number of gamma-ray bursts observed in N different flux bins. Before collecting the 
data, we think of x as a random variable with some probability distribution L(x; 9), which 
depends in some known way on a vector of M model parameters = (61,62, 0m)- 

Such model parameters might for instance be the spectral index of density fluctuations, 
the Hubble constant h, the cosmic density parameter ft or the mean redshift of gamma-ray 
bursts. We will let 9q denote the true parameter values and let 9 refer to our estimate of 
9. Since 9 is some function of the data vector x, it too is a random variable. For it to be a 
good estimate, we would of course like it to be unbiased, i.e., 

(9) = 0o, (41) 
and give as small error bars as possible, i.e., minimize the standard deviations 

Ae a ^((6l)-(6 a )Y 2 - (42) 
In statistics jargon, we want the BUE 6 a , which stands for the "Best Unbiased Estimator". 
A key quantity in this context is the so-called Fisher information matrix, defined as 

/ d 2 C \ 

F ^W^/' (43) 

where 

C=-hxL. (44) 

Another key quantity is the maximum likelihood estimator, or ML-estimator for brevity, 
defined as the parameter vector 0ml ^ na ^ maximizes the likelihood function L(x; 9). 

Using this notation, a number of powerful theorems have been proven (see e.g. |13 | .|12 | ): 

1. For any unbiased estimator, A6 a > l/%/F QQ (the Cramer-Rao inequality). 

2. If an unbiased estimator attaining ("saturating" ) the Cramer-Rao bound exists, it is 
the ML estimator (or a function thereof). 

3. The ML-estimator is asymptotically BUE. 

The first of these theorems thus places a firm lower limit on the error bars that one can 
attain, regardless of which method one is using to estimate the parameters from the data. 
You won't do better, but you might do worse. 

The normal case is that the other parameters are estimated from the data as well, in 
which case, as we have seen, the minimum standard deviation rises to 

A6 a > (F- 1 )^ 2 . (45) 

This is called the marginal error, and I reemphasise that this is normally the relevant error 
to quote. 

The second theorem shows that maximum-likelihood (ML) estimates have quite a special 
status: if there is a best method, then the ML-method is the one. Finally, the third result 
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basically tells us that in the limit of a very large data set, the ML-estimate for all practical 
purposes is the best estimate, the one that for which the Cramer-Rao inequality becomes 
an equalitjj^J It is these nice properties that have made ML-estimators so popular. 

Note that conditional and marginal errors coincide if F is diagonal. If it is not, then the 
estimates of the parameters are correlated (even if the parameters themselves are uncorre- 
cted), e.g. in the example shown in Fig. 1, estimates of the baryon density parameter Qb 
and the dark energy equation of state w = p/pc 2 are strongly correlated with WMAP data 
alone, so we will tend to overestimate both f2{, and w, or underestimate both. However, the 
value of S7b in the Universe has nothing obvious to do with w - they are independent. 




Fig. 1. The expected error ellipses for cosmological parameters (as, baryon density parameter Qh, 
and dark energy equation of state w = p/ pc 2 ) from a 3D weak lensing survey of 1000 square degrees, 
with a median redshift of 1 and a photometric redshift error of 0.15. Probabilities are marginalised 
over all other parameters, except that n = 1 and a flat Universe are assumed. Dark ellipses represent 
a prior from WMAP, pale represents the 3D lensing survey alone, and the central ellipses show the 
combination (from Kitching T., priv. comm.). 



3.1 Cramer- Rao inequality: proof for simple case 

This discussion follows closely Andrew Hamilton's lectures to the Valencia Summer School 
in 2005 0. 

At the root of the Cramer-Rao inequality is the Schwarz inequality. If we have estimators 
for two parameters 6\ and 62, then it is clear that the expectation value of 

^A§ 1 + XA9 2 y^ > (46) 

where A9 a = 8 a — 0Q a is the difference between the estimate and the true value. Evidently 
equation [46] holds for any A. It is easy to show that the left- hand-side is a minimum if we 
choose A = — (A9iA02)/((A02) 2 ), from which we obtain the Schwarz inequality: 

(M)^ ((^2)^ > (A9, Ad 2 ? (47) 
3 This is sometimes called 'saturating the Cramer-Rao bound' 



Statistical techniques in cosmology 11 

Let us treat the simple case of one parameter, and data x. An unbiased estimator 8 of a 
parameter whose true value is 9q is the value 9 which satisfies 



'o) = (0- 9 )L(x;6)dx = 0. 



Differentiating with respect to 9 we get 



QT 

{9-9 )—dx + I L{x-9)dx = Q. 



The last integral is unity, and we can therefore write 



(9-9 )- 



89 



and the Schwarz inequality gives 



(48) 



(49) 



(50) 



> 



( dlnL V 

l oe ) 



(51) 



The final expression is obtained by differentiating J L(x;9)dx — 1 twice with respect to 9 
to show 

d 2 lnL /<91ni x2 
+ 



„ 8 fd\uL Ti f 



89 J 89 
so we obtain the Cramer-Rao inequality 



89 2 



89 



L(x; 9)dx 



{(0-0o?)> 



1 



/ d 2 In L \ 
\ 89 2 I 



"60 



(52) 



(53) 



Note that for a single variable the conditional error is the same as the marginal error - the 
Fisher 'matrix' has rank 1. 



Combining experiments 

If the experiments are independent, you can simply add the Fisher matrices (why?). Note 
that the marginal error ellipses (marginalising over all but two variables) in the combined 
dataset can be much smaller than you might expect, given the marginal error ellipses for 
the individual experiments, because the operations of adding the experimental data and 
marginalising do not commute. 



3.2 The Gaussian Case 



Let us now explicitly compute the Fisher information matrix for the case when the probabil- 
ity distribution is Gaussian, i.e., where (dropping an irrelevant additive constant iVm[27r]) 

2C = lndctC + (x-/x)C- 1 (x-/x) T , (54) 

where in general both the mean vector fj, and the covariance matrix 

C=((x-/,)(x-/x) T ) (55) 

depend on the model parameters 6. Although vastly simpler than the most general situation, 
the Gaussian case is nonetheless general enough to be applicable to a wide variety of problems 
in cosmology. Defining the data matrix 

D = (x-/x)(x-/x) T (56) 
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and using the matrix identity (see exercises) lndet C = TrlnC, where Tr indicates trace, we 



can re- write ( 54 ) as 

2£ = Tr[lnC + C- 1 D] . (57) 
We will use the standard comma notation for derivatives, where for instance 

C '«4 C (58) 

Since C is a symmetric matrix for all values of the parameters, it is easy to see that all the 
derivatives C, Q , C, a p, will also be symmetric matrices. Using the matrix identities (C^ 1 )^ = 
— C _1 C, Q C _1 and (lnC), a = C _1 C, Q (see exercises), we find 

2£, Q = Tr[C- 1 C, Q -(r 1 C, Q C- 1 D + C- 1 D, a ] . (59) 

When evaluating C and fj, at the true parameter values, we have (x) = fi and (xx T ) 
C + [ifi T , which gives 

' (D) = C, 




0, 

T T 



(60) 



Using this and equation ( 59 ), we obtain (C, a ) = 0. In other words, the ML-estimate is correct 



on average in the sense that the average slope of the likelihood function is zero at the point 



corresponding to the true parameter values. Applying the chain rule to equation (59), we 
obtain 



2£, a /3 — Tr[ — C C, a C 1 C, ( g+C C, 



a0 



c- x (c, a c -1 ^ +c,/g c- x c, a )c- : d 

C- 1 (C >a C- 1 D^ + C )/ ,C- 1 D ia ) 

C^C^C^D + C^D,^]. (61) 



Substituting this and equation (601 into equation (43) and using the trace identity Tr[AB] = 



Tr[BA], many terms drop out and the Fisher information matrix reduces to simply 

F Q /J = (Aa/3 ) = ^[C^C^C^C^ + C^M^], (62) 

where we have defined the matrix M a p = (D, a p ) = /j,, a /ix,J. 
The Fisher matrix requires no data. 

This result is extremely powerful. If the data have a (multivariate) gaussian distribution 
(and the errors can be correlated; C need not be diagonal), and you know how the means fi 
and the covariance matrix C depend on the parameters, you can calculate the Fisher Matrix 
before you do the experiment. The Fisher Matrix gives you the expected errors, so you know 
how well you can expect to do if you do a particular experiment, and you can then design 
an experiment to give you, for example, the best (marginal) error on the parameter you are 
most interested in. 

Note that if the prior is not uniform, then you can simply add a 'prior matrix' to the 
Fisher matrix before inversion. Fig. [T] shows an example, where a prior from CMB experi- 
mental results has been added to a hypothetical 3D weak lensing survey. 

Treat the Fisher errors as a one-way test: you might not achieve errors which are this 
small, but you won't do better. So if you want to measure some quantity with an accuracy 
of a metre, and a Fisher analysis tells you the error bar is the size of Belgium, give up. 

Finally, note that this analysis assumes that the data do not depend on the parameters. 
Normally this is the case - you simply measure the data, right? Be careful - if you construct a 
new dataset (e.g. estimates of a galaxy correlation function), you may need to assume what 
the parameters are (e.g. what the cosmology is, to calculate the pair separations). In this 
case, you should assume a fiducial set of cosmological parameters to calculate your dataset, 
and not change the data as you change the parameters. You then need to do more work on 
the theory side, but the error bars should be robust. 



Reparametrisation 
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Finally we mention that sometimes the gaussian approximation for the likelihood surface 
is not a very good approximation. With a good theoretical model, it is possible to make 
nonlinear transformations of the parameters to make the likelihood more gaussian. See |14j 
for more details. 

iCosmo: a great resource 

icosmo.org has a web-based calculator for Fisher matrices for cosmology. You can download 
results, or even the source code, to compute Fisher matrices for various experiments in 
lensing, BAOs, supernovae etc, or custom-design your own survey. It also gives many other 
useful things, such as lensing power spectra, angular diameter distances etc. 



4 Numerical methods 

If the problem has only two or three parameters, then it may be possible to evaluate the 
likelihood on a sufficiently fine grid to be able to locate the peak and estimate the errors. If 
the dimensionality of the parameter space is very large, then, as the number of grid points 
grows exponentially with dimension, it becomes rapidly unfeasible to do it this way. In fact, 
it's very inefficient to do this anyway, as typically most of the hypervolume has very small 
likelihood so is of little interest. There are various ways to sample the likelihood surface 
more efficiently, concentrating the points more densely where the likelihood is high. We 
cover here the most common method (MCMC), and a relatively new method (to cosmology), 
Hamiltonian Monte Carlo, which could take over, as it seems more efficient in cases studied. 

4.1 Monte Carlo Markov Chain (MCMC) method 

The aim of MCMC is to generate a set of points in the parameter space whose distribution 
function is the same as the target density, in this case the likelihood, or more generally 
the posterior. MCMC makes random drawings, by moving in parameter space in a Markov 
process - i.e. the next sample depends on the present one, but not on previous ones. By 
design, the resulting Markov Chain of points samples the posterior, such that the density of 
points is proportional to the target density (at least asymptotically), so we can estimate all 
the usual quantities of interest from it (mean, variance, etc). The number of points required 
to get good estimates is said to scale linearly with the number of parameters, so very quickly 
becomes much faster than grids as the dimensionality increases. In cosmology, we are often 
dealing with around 10-20 parameters, so MCMC has been found to be a very effective tool. 

The target density is approximated by a set of delta functions (you may need to nor- 
malise) 



The basic procedure to make the chain is to generate a new point 9* from the present 
point 8 (by taking some sort of step), and accepting it as a new point in the chain with a 
probability which depends on the ratio of the new and old target densities. The distribution 
of steps is called the proposal distribution. The most popular algorithm is the Metropolis- 
Hastings algorithm, where the probability of acceptance is 




(63) 



from which we can estimate any integrals (such as the mean, variance etc.): 




(64) 
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p(acceptance) 



x p{o*) q {e*\e) 



(65) 



p(O)q(8\0*) 

where the proposal distribution function is q{6 * \ ff) for a move from 6 to 6* . 

If the proposal distribution is symmetric (as is often the case) , the algorithm simplifies 
to the Metropolis algorithm: 

p(6*) 



p(acceptance) = min 



1; p{9) 



(66) 



• Choose a random initial starting point in parameter space, and compute the target 
density. 

• Repeat: 

• Generate a step in parameter space from a proposal distribution, generating a new trial 
point for the chain. 

• Compute the target density at the new point, and accept it (or not) with the Metropolis- 
Hastings algorithm. 

• If the point is not accepted, the previous point is repeated in the chaii^ 

• End Repeat: 

The easy bits: this is trivial to code - you might just take a top-hat proposal distribution 
in each parameter direction, and it should work. The harder parts are (even in the tophat 
case): choosing an efficient proposal distribution; dealing with burn-in and convergence. 



Proposal distribution 



If the proposal distribution is small, in the sense that the typical jump is small, then the 
chain may take a very long time to explore the target distribution, and it will be very 
inefficient. Since the target density hardly changes, almost all points are accepted, but it 
still takes forever. This is an example of poor mixing. If the proposal distribution is too 
large, on the other hand, then the parameter space is explored, but the trial points are often 
a long way from the peak, at places where the target density is low. This is very inefficient as 
well, since they are almost always rejected by the Metropolis-Hastings algorithm. So what is 
best? You might expect the chain to do well if the proposal distribution is 'about the same 
size as the peak in the target density', and you would be right. In fact, a very good option 
is to draw from a multivariate gaussian with the Fisher matrix as Hessian. However, having 
said that, if you want something quick to code (but sub-optimal), a top- hat of about the 
right dimensions will do a decent job. You might have to run a preliminary chain or two to 
get an idea of the size of the target, but the rules say you are not allowed to change the 
proposal distribution in a chain, so once you have decided you must throw away your chains 
and begin again. 



Burn-in and convergence 



Theory indicates that the chain should fairly sample the target distribution once it has 
converged to a stationary distribution. This means that the early part of the chain (the 
'burn-in' are ignored, and the dependence on the starting point is lost. It is vitally important 
to have a convergence test. Be warned that the points in a MCMC chain are correlated, and 
the chain can appear to have converged when it has not (one can reduce this problem by 
evaluating the correlation function of the points, and 'thinning' them by (say) taking only 
every third point (or whatever is suggested by the correlation analysis) . Fig. [2] shows one 
such example. The classic test is the Gelman-Rubin (1992) convergence criterion. Start M 



4 It is a common mistake to neglect to do this 
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chains, each with 27V points, starting at well-separated parts of parameter space. In this 
example, the first N are discarded as burn-in. 

The idea is that you have two ways to estimate the mean of the parameters - either treat 
the combined chains as a single dataset, or look at the means of each chain. If the chains 
have converged, these should agree within some tolerance. 

Following [24] . let Oj represent the point in parameter space in position / of chain J. 
Compute the mean of each chain (J): 



JV 

F = 

N 



i=i 



and the mean of all the chains 



N M 



i=i j=i 

The chain-to-chain variance B is 

M 

(M — 1) j=i 
and the average variance of each chain is 

N M 

M(N 



n =^ £ny E£(# (70) 



i=i j=i 

Under convergence, W and B/N should agree. 
The weighted estimate of the variance, 

9 N- 1 B 

overestimates the true variance if the starting distribution is overdispersed. 
Accounting for the variance of the means gives an estimator of the variance 

V is an overestimate, and W is an underestimate. The ratio of the two estimates is 

a = m±M±i). (73) 

R should approach unity as convergence is achieved. How close to R = 1? Opinions differ; 
I have seen suggestions to run the chain until the values of R are always < 1.03, or < 1.2, 
but a proof would be nice. 



CosmoMC 

Excellent resource. MCMC sampler with cosmological data (CMB + support for LSS, SNe). 
|http: / / cosmologist.info/cosmomc/| 

For more details on MCMC in general, see [51 UM |2"5] . 



4.2 Hamiltonian Monte Carlo 



As we've seen, the proposal distribution has to be reasonably finely tuned to ensure good 
mixing, but not be too inefficient. What we would really like is to be able to take rather big 
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Fig. 2. Examples of unconverged chains. The left panel suggests the chain has converged, but 
the right panel shows the chain and a second one, also apparently converged, but showing clearly 
different distributions. From Verde et al. ApJS, 148, 195 (2003). 



jumps, so the chain is well-mixed, but in such a way that the probability of acceptance of 
each point is still high. Hamiltonian (or Hybrid) Monte Carlo (HMC) tries to do this, via a 
rather clever trick. In practice, it seems that it can be about M times faster than MCMC 
in M dimensions. Typical applications report that 4 times shorter chains give the same 
accuracy as MCMC. Originally developed for particle physics [3], there is a nice exposition 
in astrophysics by Hajian [7]. 

HMC works by sampling from a larger parameter space than we want to explore, by 
introducing M auxiliary variables, one for each parameter in the model. To see how this 
works, imagine each of the parameters in the problem as a coordinate. HMC regards the 
target distribution which we seek as an effective potential in this coordinate system, and 
for each coordinate it generates a generalised momentum, i.e. it expands the parameter 
space from its original 'coordinate space' to a phase space, in which there is a so-called 
'extended' target distribution. It tries to sample from the extended target distribution in 
the 2M dimensions. It explores this space by treating the problem as a dynamical system, 
and evolving the phase space coordinates by solving the dynamical equations. Finally, it 
ignores the momenta (marginalising, as in MCMC), and this gives a sample of the original 
target distribution. 

There are several advantages to this approach. One is that the probability of acceptance 
of a new point is high - close to unity; secondly, the system can make big jumps, so the 
mixing is better and the convergence faster. In doing the big jumps, it does some addi- 
tional calculations, but these do not involve computing the likelihood, which is typically 
computationally expensive. 

Let us see how it works. If the target density in M dimensions is p(9), then we define a 
potential 

£7(0) = -\np(0). (74) 

For each coordinate 9 a , we generate a momentum u a , conveniently from a normal distribu- 
tion with zero mean and unit variance, so the M-dimensional momentum distribution is a 
simple multivariate gaussian which we denote A/"(u) . We define the kinetic energy 



K(U) EE lu T U, 



(75) 
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and the Hamilitonian is 

H(0,u) = U(0) + K{u). (76) 
The trick is that we generate chains to sample the extended target density 

p{0,u)=exp[-H(0,u)]. (77) 

Since this is separable, 

p(6, u) = exp {-U{d)} exp [-K{n)\ oc p(0)JV(u) (78) 

and if we then marginalise over u by simply ignoring the u coordinates attached to each 
point in the chain (just as in MCMC), the resulting marginal distribution samples the desired 
target distribution p{0). This is really very neat. 

The key is that if we exactly solve the Hamiltonian equations 

*■ - - *u < 79 > 

then H remains invariant, so the extended target density is always the same, and the ac- 
ceptance is unity. Furthermore, we can integrate the equations for a long time if we wish, 
decorrelating the points in the chain. 

There are several issues to consider. 

• We seem to need the target density to define the potential, but this is what we are looking 
for. We need to approximate it. 

• The aim is to do this fast, so we do not want to do many operations before generating a 
new sample. An easy way to achieve this is to employ a simple integrator (e.g. leap-frog; 
even Euler's method might do) and take several fairly big steps before generating a new 
point in the chain. 

• We need to ensure we explore the extended space carefully. 

The result of the two approximations is that H will not be quite constant. We deal with 
this as in MCMC by using the Metropolis algorithm, accepting the new point (9* , u*) with 
a probability 

min{l,exp[-H(0*,u*) + H(0,u)]}, (80) 
otherwise we repeat the old point (0, u) as usual. 

How do we approximate? We want the gradients to be cheap to compute. It is usual to 
generate an approximate analytical potential by running a relatively short MCMC chain, 
computing the covariance of the points in the chain, and approximating the distribution by 
a multivariate gaussian by a multivariate gaussian with the same covariance. The gradients 
are then consequently easy to compute analytically. 

The last point is that if we change the momentum only with Hamilton's equations of 
motion, we will restrict ourselves to a locus in phase space, and the target distribution will 
not be properly explored. To avoid this, a new momentum is generated randomly when each 
point in the chain is generated. The art is to choose a good step in the integration, and 
the number of steps to take before generating a new point. Perhaps unsurprisingly, choosing 
these such that the new point differs from the previous one by about the size of the target 
peak works well. Thinning can be performed, and a convergence test must still be applied. 
Fig. [3] shows a comparison between HMC and MCMC for a simple case of a 6D gaussian 
target distribution, from [7]. 



5 Model Selection 



Model selection is in a sense a higher-level question than parameter estimation. In parameter 
estimation, one assumes a theoretical model within which one interprets the data, whereas 
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Fig. 3. A comparison of HMC sampling (top) and MCMC sampling (bottom). Note that the 
computer time required to generate each point in the HMC sampling will be larger than that of 
MCMC, so the actual gains are less than appears. From [7]. 



in model selection, one wants to know which theoretical framework is preferred, given the 
data (regardless of the parameter values). The models may be completely different (e.g. 
compare Big Bang with Steady State, to use an old example), or variants of the same 
idea. E.g. comparing a simple cosmological model where the Universe is assumed flat and 
the perturbations are strictly scale-invariant (n — 1), with a more general model where 
curvature is allowed to vary and the spectrum is allowed to deviate from scale-invariance. 
The sort of question asked here is essentially 'Do the data require a more complex model?'. 
Clearly in the latter type of comparison \ 2 itself will be of no use - it will always reduce if 
we allow more freedom. There are frequentist ways to try and answer these questions, but 
we are all by now confirmed Bayesiansj^J so will approach it this way. 



5.1 Bayesian evidence 

The Bayesian method to select between models (e.g. General Relativity, & modified gravity) 
is to consider the Bayesian evidence ratio. Essentially we want to know if, given the data, 
there is evidence that we need to expand the space of gravity models beyond GR. Assuming 
non-commital priors for the models (i.e. the same a priori probability), the probability of 
the models given the data is simply proportional to the evidence. 

We denote two competing models by M and M' . We assume that M' is a simpler model, 
which has fewer (n' < n) parameters in it. We further assume that it is nested in Model M', 
i.e. the n' parameters of model M' are common to M, which has p = n—n' extra parameters 
in it. These parameters are fixed to fiducial values in M' . 

We denote by x the data vector, and by 6 and 6' the parameter vectors (of length n and 

n>). 

Apply Rule 1: Write down what you want to know. Here it is p(M|x) - the probability 
of the model, given the data. 

The posterior probability of each model comes from Bayes' theorem: 

p(x\M)p(M) 

p(M|x) = — — ( } 

and similarly for M' . By marginalisation p(x|M ), known as the Evidence, is 

p{x\M) = J d0p(x\0M)p(0\M), (82) 

which should be interpreted as a multidimensional integration. Hence the posterior relative 
probabilities of the two models, regardless of what their parameters ar^J is 

5 If not, please leave the room 

6 If a model has no parameters, then the integral is simply replaced by p(xjM) 
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p(M'|x) = p(M>) fdO'p(x\9'M')p(0'\M>) 

p(M|x) p(M) Jd9p(x\9M)p(9\M) ' 1 j 

With non-committal priors on the models, p(M') = p[M), this ratio simplifies to the ratio 
of evidences, called the Bayes Factor, 

= 5de'p{x.\9'M>)p{e'\M>) 

~ J d6p{x\6M)p[0\M) ' 1 ' 

Note that the a complicated model M will (if M' is nested) inevitably lead to a higher 
likelihood (or at least as high), but the evidence will favour the simpler model if the fit is 
nearly as good, through the smaller prior volume. 

We assume uniform (and hence separable) priors in each parameter, over ranges A9 (or 
AO'). Hence p{6\M) = (A9 1 . . . A9 n )~ 1 and 

_ Jd8'p(x\9',Mi) A9 1 ...A9 n 

fdOp(x\0,M) A9 , x ...A9 l n ,' 1 ' 

Note that if the prior ranges are not large enough to contain essentially all the likelihood, 
then the position of the boundaries would influence the Bayes factor. In what follows, we 
will assume the prior range is large enough to encompass all the likelihood. 

In the nested case, the ratio of prior hypervolumes simplifies to 

= A6 n , +1 . . . A9 n > +P , (86) 

where p = n — n! is the number of extra parameters in the more complicated model. 

Here we see the problem. The evidence requires a multidimensional integration over the 
likelihood and prior, and this may be very expensive to compute. There are various ways 
to simplify this. One is analytic - follow the Fisher approach and assume the likelihood is 
a multivariate gaussian, others are numerical, such as nested sampling, where one tries to 
sample the likelihood in an efficient way. There are others, but we will focus on these. Note 
that shortcuts with names such as AIC and BIC may be unreliable as they are based on the 
best-fit x 2 , and from a Bayesian perspective we want to know how much parameter space 
would give the data with high probability. See |17) for more discussion. 



5.2 Laplace approximation 



The Bayes factor in equation ( 85 ) still depends on the specific dataset x. For future experi- 



ments, we do not yet have the data, so we compute the expectation value of the Bayes factor, 
given the statistical properties of x. The expectation is computed over the distribution of x 
for the correct model (assumed here to be M). To do this, we make two further approxima- 
tions: first we note that B is a ratio, and we approximate (B) by the ratio of the expected 
values, rather than the expectation value of the ratio. This should be a good approximation 
if the evidences are sharply peaked. 

We also make the Laplace approximation, that the expected likelihoods are given by 
multivariate Gaussians. For example, 



(p(x|0,M)) =L exp 



1 



(9 — 9 ) a V a p{9 — 9 )p 



(87) 



and similarly for (p(x|0 , M')). This assumes that a Taylor expansion of the likelihood around 
the peak value to second order can be extended throughout the parameter space. F a p is the 



Fisher matrix, given for Gaussian-distributed data by equation (62): 



1 



'a/3 



Tr [C 1 C ) „C 1 C^ + C 1 (fi,pfJ? a +M,a^)] 
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C is the covariance matrix of the data, and /i its mean (no noise) . Commas indicate partial 
derivatives with respect to the parameters. For the correct model M, the peak of the expected 
likelihood is located at the true parameters 9q. Note, however, that for the incorrect model 
M', the peak of the expected likelihood is not in general at the true parameters (see Fig. [4] 
for an illustration of this). This arises because the likelihood in the numerator of equation 



(85) is the probability of the dataset x given incorrect model assumptions. 



The Laplace approximation is routinely used in forecasting marginal errors in parameters, 
using the Fisher matrix. Clearly the approximation may break down in some cases, but for 
Planck, the Fisher matrix errors are reasonably close to (within 30% of) those computed 
with Monte Carlo Markov Chains. 

If we assume that the posterior probability densities are small at the boundaries of the 
prior volume, then we can extend the integrations to infinity, and the integration over the 
multivariate Gaussians can be easily done. This gives, for M, (2 7 r) n / 2 (detF)- 1 / 2 , so for 
nested models, 

(B) = (2n)-^^^^A9 n , +1 . . . A9 n , +P . (89) 

An equivalent expression was obtained, using again the Laplace approximation by |15j . The 
point here is that with the Laplace approximation, one can compute the L' /Lq ratio from 
the Fisher matrix. To compute this ratio of likelihoods, we need to take into account the fact 
that, if the true underlying model is M, in M' (the incorrect model), the maximum of the 
expected likelihood will not in general be at the correct values of the parameters (see Fig. |4| . 
The n' parameters shift from their true values to compensate for the fact that, effectively, 
the p additional parameters are being kept fixed at incorrect fiducial values. If in M', the 
additional p parameters are assumed to be fixed at fiducial values which differ by Sip a from 
their true values, the others are shifted on average by an amount which is readily computed 
under the assumption of the multivariate Gaussian likelihood: 

59' a = -(F'- 1 ) a pGp C 6i) ( a,P = l...ri,C = l...p (90) 

where 

G K = ^Tr [C-^pC- 1 ^ + C" 1 ^,^ + fi,^)] , (91) 

which we recognise as a subset of the Fisher matrix. For clarity, we have given the additional 
parameters the symbol ip^; £ = l...p to distinguish them from the parameters in M'. 




Fig. 4. Illustrating how assumption of a wrong parameter value can influence the best-fitting value 
of other model parameters. Ellipses represent iso-likelihood surfaces, and here in the simpler model, 
the parameter on the horizontal axis is assumed to take the value given by the vertical line. Filled 
circles show the true parameters in the more complicated model, and the best-fit parameters in the 
simpler model. From [5]. 



With these offsets in the maximum likelihood parameters in model M' , the ratio of 
likelihoods is given by 
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L'n Uvx V [ ~6e a F afj S0p) (92) 



where the offsets are given by S9 a = S0' a for a < n' (equation 90), and 80 a = Sip a ^ n i for 
a > n' . 

The final expression for the expected Bayes factor is then 

(B) = ^y P/2 ^J== ex P (-\sO a ? aP SB^) f[ A9 n , +q . (93) 



q=l 



Note that F and F _1 are n x n matrices, F' is n' X n', and G is an nf x p block of the full 
n x n Fisher matrix F. The expression we find is a specific example of the Savage-Dickey 
ratio (see e.g. [13]). For a nested model with a single additional parameter # i5 

(B) = «f (94) 

Here we explicitly use the Laplace approximation to compute the offsets in the parameter 
estimates which accompany the wrong choice of model, and compute the evidence ratio 
explicitly. Finally, note that this is the expected evidence ratio (nearly); it does not address 
the issue of what the distribution of evidence ratios should be. 

Note that the 'Occam's razor' term, common in evidence calculations, is to some extent 
encapsulated in the term (27r) _p / 2 ^=|=, multiplied by the prior product: models with more 
parameters are penalised in favour of simpler models, unless the data demand otherwise. In 
cases where the Laplace approximation is not a good one, other techniques must be used, 
at more computational expense. 

It is perhaps worth remarking that Occam's razor appears not to be fully incorporated 
into this term, as can be seen by considering a situation where the data do not depend 
at all on an additional parameter. In this case, the Bayesian evidence ratio is unity, so 
disappointingly no preference is shown at all. 

As an example, consider testing General Relativity against other gravity theories, which 
predict a different growth rate of perturbations, d\n5/d\n fl m — 7, where 7 = 0.55 for GR, 
and (for example) 7 = 0.68 for a flat DGP braneworld model. This can be probed with weak 
lensing, and we ask the question do the data favour a model where 7 is a free parameter, 
rather than being fixed at 0.55? 

We take a prior range A^ = 1, and we ask the question of how different the growth 
rate of a modified-gravity model would have to be for these experiments to be expected to 
favour a relaxation of the gravity model from General Relativity. This is shown in Fig(5] It 
shows how the expected evidence ratio changes with progressively greater differences from 
the General Relativistic growth rate. We see that a next-generation weak lensing survey 
could even distinguish 'strongly' £7 = 0.048. Note that changing the prior range A'j by a 
factor 10 changes the numbers by ~ 0.012, so the dependence on the prior range is rather 
small. 

If one prefers to ask a frequentist question, then a combination of WL+PZaracfc+BAO+SN 
should be able to distinguish £7 = 0.13, at 10.6cr. Alternatively, one can calculate the 
expected error on 7 [T] within the extended model M. In this section, we are asking a 
slightly different question of whether the data demand that a wider class of models needs 
to be considered at all, rather than estimating a parameter within that wider class. 

The case for a large, space-based 3D weak lensing survey is strengthened, as it offers the 
possibility of conclusively distinguishing Dark Energy from at least some modified gravity 
models. 



5.3 Numerical sampling methods 



In order to compute the evidence numerically, the prior volume needs to be sampled, in 
much the same way as in parameter estimation, except the requirements are slightly more 
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Fig. 5. The expected value of | ln(B)| from a future large-scale deep weak lensing survey, as might 
be done with Euclid, JDEM or LSST, in combination with CMB constraints from Planck, as a 
function of the difference in the growth rate between the modified-gravity model and General 
Relativity. The crossover at small 5"/ occurs because Occam's razor will favour the simpler (General 
Relativity) model unless the data demand otherwise. To the left of the cusp, GR would be likely to 
be preferred by the data. The dotted vertical line shows the offset of the growth factor for the flat 
DGP model. The descriptors are the terminology of Jeffreys (1961) [11]. From [9] 

stringent. MCMC could be used, although there are claims that it is not good at exploring 
the parameter space adequately. I find these claims puzzling, as in evidence calculations 
we are doing an TV-dimensional integral, which is not so different from doing the (N — 1)- 
dimensional integral in MCMC to get marginal errors. However, here are a couple of other 
sampling techniques. 

The VEGAS algorithm, with rotation 

This is suitable for single-peaked likelihoods. It is in Numerical Recipes, but needs a modifi- 
cation for efficiency. Essentially, one seeks to sample from a distribution which is close to the 
target distribution (sampling from a different distribution is called importance sampling), 
but one does not know what it is yet. One can do this iteratively, sampling from the prior 
first, then estimating the posterior to get a first guess at the posterior, and using that to 
refine the sampling distribution. 

Now, one does not want to draw randomly from a non-separable function of the parame- 
ters, as a moment's thought will tell you that this is computationally very expensive, so one 
seeks a separable function, so one can then draw the individual parameters one after the 
other from N distributions. This works well if the target distribution is indeed separable in 
the parameters, but not otherwise. 

The key extra step |20) is to rotate the parameter axes, which can be done by computing 
(for example) the moments of the distribution (essentially the moment-of-inertia) after any 
step, and diagonalising it to find the eigenvectors. 

The probability to be sampled is 



p(6)ccg 1 (6 1 )g 2 (0 2 )...g M (e M ). 



(95) 



where it can be shown that 




(96) 



where / is the desired target distribution. Note that g depends on all other gs. g can be 
improved iteratively. 
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DENSITY OF SAMPLING POINTS 




0.1 D.15 0.2 0.25 0.3 0.35 



Fig. 6. The VEGAS Sampling algorithm, applied to supernova data (from [20] 
Nested sampling 

Other sampling methods can also be very effective (see e.g. [23l O [lOj [18] ) . Nested sampling 
was introduced by Skilling [5T]. One samples from the prior volume, and gradually concen- 
trates more points near the peak of the likelihood distribution, by repeatedly replacing the 
point with the lowest target density by one drawn strictly from the prior volume with higher 
target density. This has proved effective for cosmological model selection. I will not go into 
details here, except to say that the key is in drawing the new point from a suitable subset 
of the prior volume. This must be increasingly smaller as the points get more confined, 
otherwise the trial points will be rejected with increasing frequency and it becomes very 
inefficient, but it must also be a large enough subset that the entire prior volume above the 
lowest target density is almost certainly included inside. For further details, see the origi- 
nal papers. Note that, for multimodal target distributions, a modification called MultiNest 
exists [4]. CosmoNest and MultiNest are publicly-available additions to CosmoMC. 



6 Summary 

We have explored the use of (principally) Bayesian methods for parameter estimation and 
model selection, and looked at some common numerical methods for performing these tasks. 
We have also shown how it is possible to work out the (minimum) errors on parameters, 
and on model probabilities, in advance of doing an experiment, to see if it is worthwhile. 
These are the common Fisher matrix approach for parameter estimation, and the expected 
evidence for model selection, which requires no more than the Fisher matrix to compute. 



7 Exercises 

1. A Bayesian Exercise. This is a variant of the famous Monty Hall problem, named after a 
game show host. For this exercise, you must formulate the problem in a fully Bayesian way. 

You are participating in a game show where there are three small doors presented to 
you. You are told that behind two of them there is a bottle of fine italian wine, and behind 
the third is a can of the famous Scottish soft drink, Irn Bru. You will win one of them, and 
naturally you prefer the Irn Bru. The game works as follows: 
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a) You point to, but do not open, one of the doors. 

b) The game show host opens one of the other doors, revealing a bottle of wine. 

c) You now either open the door which you originally selected, or switch to the third door. 
Which should you open? 

2. This is a model selection problem, and will get you thinking about suitable priors. You 
draw lottery ticket number 1475567, and it wins. What is the probability that the lottery 
was rigged, so that the winning ticket was predetermined? 

3. An exercise in linear algebra 

a) Prove that (C -1 )^ = -C _1 C, Q C" 1 

b) Prove that (lnC), Q = C _1 C, a . 

c) Prove that In det C = Tr In C. 

Hints: 

a) What is CC" 1 ? 

b) For matrices, exp(A) = I + A + A 2 /2! + .... In is the inverse of exp. 

c) If you rotate coordinate axes (in parameter space inFcvidcncc this case), the determinant 
and trace don't change. 

4. A survey is proposed, to determine the mean number density of a certain type of astro- 
nomical object, whose positions are random. The survey measures the number of objects in 
each of N independent cells of the same size and shape, such that the mean number per cell, 
n, is 3> 1. If the cells are observed to have n, objects in them [i = 1, . . . N), show that the 
Fisher matrix (a scalar in this case) for n is 

_ N N 

~ m 

Where is most of the information coming from, the dependence of fi on n, or C? 

5. Binomial drinking - an interesting paradox. Precursor question. If we have an experiment 
with some discrete outcomes n = 1, . . . , oo, each occurring with probability P n . Argue that 
the probability of the sum of two drawings z = m + n from the distribution is 



(97) 



-1 



Every minute, as the second hand on a large clock on the wall reaches the top, I think 
about drinking a bottle of Irn Bru. I drink one with a probability p. Show that the probability 
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of the next drink being taken at the M th opportunity is 



P M = PQ™- 1 



where q = 1 — p. 



Show that the expectation value of the number of time steps between outbursts is M = 
1/p. (Hint: expand (1 — q)~ 2 in a Taylor expansion). 

Now, my friend comes in at a random time (not necessarily as the second hand reaches 
the top of the clock). Argue that the probability that I took my last drink m time steps 
previously was 

P' m =pq m - 1 - 

Show that the probability for the variable S = M + m (S > 2) is 

M=l 

By expanding a different power of 1 — q, show that the expectation value of S is 

s = 2 , 
p 

and hence that the average time between the last output and the next one is 

P 

For p = 0.1, t = 19, compared with 10 for the mean time between drinks. How is the paradox 
resolved? 
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8 Solutions to selected exercises 



1. The aim here is to do a Bayesian calculation, not necessarily to find the easiest or best 
solution. Let the doors be labelled a, b, c, where a is the door you choose initially, and b is 
the door which is opened. Many of the probabilities below should be interpreted as 'given 
that you have chosen a\ but we won't write this explicitly. 

Let p(a) = probability that a leads to the Irn Bru. p(b), p(c) similarly. 

Let B be the event that door b gets opened and leads to a bottle of wine. 

What you want is the probability that a leads to the Irn Bru, given that b is opened and 
leads to a bottle of wine. i.e. the aim is to calculate 

p(a\B). 

We can use Bayes' theorem for this: 

p(a,B) _p(B\a)p(a) 



p{a\B) 



p{B) p(B) 



Now, clearly p(a) = p(b) = p(c) = 1/3 (all doors are equally likely, before any experiment is 
done). 

p(B\a) = probability that door b is opened, given that a leads to the Irn Bru. 

Since the host could have opened either door b or c, since they both lead to wine, we 
have 

P(B\a) = \ : 

What about p{B)l It is the sum of all the joint probabilities: 

p(B) = p(B, a) + p(B, b) + p(B, c) = p(B\a)p(a) + p(B\b)p(b) + p(B\c)p(c), 

each of which we can calculate. p(a) = p(b) = p(c) = 1/3, as before, and p(B\a) = 1/2 as 
before. Now 

p(B\b) = : 

The host will not open b since it leads to the Irn Bru in this case. 

p(B\c) is the most interesting. Given that you have chosen a (remember this is implicit 
throughout), then if c leads to the Irn Bru, then the host must open door b, i.e. 

p(B\c) = 1 

So the probability that your original choice a leads to the Irn Bru is 

p(B\a)p(a) 

P{alB > = p(B\a)p(a)+p(B\b)p(b)+p(B\c)p(c) (98) 
l l 

_ 2 3 



" l§ + (0x§) + (lx|) 
1 

~ 3 

So you would double your chances (from 1/3 to 2/3) if you switch to the other door. 

PS There is a really easy way to do this by lateral thinking. After the host has opened a 
door, there is one remaining door with Irn Bru, and one door with wine. Hence if you swap, 
you will either change success into failure, or failure into success. Since your original chances 
of success were only 1/3, you improve this to 2/3 if you swap. Easy. 

3. (a) Prove that (C^ 1 )^ = -C _1 C, a C" 1 

Since CC _1 = I, its derivative is zero. Hence C(C _1 ) iCt + C iCt C _1 = 0. Result follows after 
premultiplication by C _1 . 
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(b) Prove that (lnC), Q = C^C, 
Let A = InC, so 



A 2 A" 
C = expA=l + A+- + ... = ^- 

n=0 



Hence 



ln-1 



A" 



ra=l V m=0 

and result follows. 

(c) Prove that lndet C = Trln C. 

C is a symmetric matrix and therefore be diagonalised. In the new basis (where we 
have a new set of parameters which are linear combinations of the old ones), the diagonal 
components of C are strictly positive (they are variances of the new parameters). 

Since the trace and determinant are unchanged by the diagonalisation, we can prove 
the result in the rotated system. If C is diagonal, then InC is diagonal, with components 
lnCn, lnC 22 , .... So Trln C = ^ n lnC m . Since detC = J]„ C nn , 

In det C = ^ In C„„ = Tr In C. 

n 

This proof is rigorous, but there may be a neater solution without diagonalisation. 

4. For n ^> 1, we can approximate the Poisson distribution by a gaussian, with (rij) = n, and 
of = n. Hence /x T = n(l, 1, ... 1) and C = ndiag(l, 1, . . . , 1). Hence Ci = diag(l, 1, . . . , 1), 
and Mij is an N x N matrix filled with 2s. The result 



follows. The first term arises from C,i, and the second from fi 1 . The second dominates since 
n > 1. 



